Test

library(raster)

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

library(neonAOP)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(rhdf5)
library(rgeos)

## rgeos version: 0.3-19, (SVN revision 524)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-3 
##  Polygon checking: TRUE

library(ggplot2)

Load canopy height model

chm <- raster("../NEONdata/D17-California/SOAP/2013/lidar/SOAP_lidarCHM.tif")

Identify plot boundaries

We wanted to match up our in situ data with the hyperspectral data so that we could subset a hyperspectral flight line. As a first step, we identified plot boundaries based on the stem locations in the in situ data.

First, we load the in situ data as a shapefile

stem.map <- readOGR("../NEONdata/D17-California/SOAP/2013/insitu/veg-structure",
                    "soap_stems")

## OGR data source with driver: ESRI Shapefile 
## Source: "../NEONdata/D17-California/SOAP/2013/insitu/veg-structure", layer: "soap_stems"
## with 1231 features
## It has 26 fields

# look at a plot
plot(chm)
plot(stem.map, add=TRUE)

Use stem locations to create plot boundaries

# group stems by plotid, record the max and min northing and easting values
# this will be used later to create a shapefile for plot boundaries
stem.map.extent <- stem.map@data %>% 
  group_by(plotid) %>%
  summarise(northing.max = max(northing) + 5,
            northing.min = min(northing) - 5,
            easting.max = max(easting) + 5,
            easting.min = min(easting) - 5)

# assign new variables for use with previously created code
yPlus <- stem.map.extent$northing.max
yMinus <- stem.map.extent$northing.min
xPlus <- stem.map.extent$easting.max
xMinus <- stem.map.extent$easting.min

# code from NEON tutorial on creating square plot extents
square <- cbind(xMinus, yPlus, 
                xPlus, yPlus, 
                xPlus, yMinus, 
                xMinus, yMinus, 
                xMinus, yPlus)

ID <- stem.map.extent$plotid

Create spatial polygons using the coordinates

# Create a function to do this
polys <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)  # take a list and create a matrix
  Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID),proj4string=CRS(as.character("+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")))

Create shapefile

polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

Plot this with our CHM

plot(chm)
plot(polys.df, add=TRUE)

Look at all the hyperspectral flightlines

We used Leah’s code to look through all the hyperspectral flightlines. I’m not going to rehash it all here, but we decided to narrow down to one flightline that covered 4 plots.

We also used Leah’s code to get extents for all flightlines and saved this on our local computer.

Set data directory to access hyperspectral flightline from hard drive

## SOAP Clip
# the name of the site
site <- "SOAP"
domain <- "D17"
fullDomain <- "D17-California"
level <- "L1"
dataType <- "Spectrometer"
level <- paste0(site,"_L1")
year <- "2013"
productType <- paste0(site,"_", dataType)
dataProduct <- "Reflectance"

drivePath <- "Volumes"
driveName <- "AOP-NEON1-4"

dataDir <- file.path(drivePath, driveName,
                      domain,
                      site, year, level, productType, dataProduct)
dataDir <- paste0("/", dataDir)

Import flightline

The right boundary of the flightline appears on the plot below

flight1 <- readOGR("exports/SOAP_flightLines","NIS1_20130612_104651_atmcor")

## OGR data source with driver: ESRI Shapefile 
## Source: "exports/SOAP_flightLines", layer: "NIS1_20130612_104651_atmcor"
## with 1 features
## It has 1 fields

# look at this with our plots
plot(chm)
plot(polys.df, add=TRUE)
plot(flight1, add=TRUE)

Identify plots within flightline

Note that 4 plots are close to the center of the flightline. We will subset for those 4 plots.

# choose the plots that intersect with flight 1 for extracting HSI
flight1.plots <- intersect(polys.df, flight1)
flight1.plots

## class       : SpatialPolygonsDataFrame 
## features    : 4 
## extent      : 297029, 297132.1, 4100078, 4100919  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       :    id.1,                           id.2 
## min values  : SOAP139, NIS1_20130612_104651_atmcor.h5 
## max values  : SOAP555, NIS1_20130612_104651_atmcor.h5

# check this subset
plot(chm)
plot(flight1.plots, add=TRUE)

Boundary for hyperspectral subset

Now we create a boundary that includes all 4 plots to subset the hyperspectral data.

# thanks for the code, leah!

# define the CRS definition by EPSG code
epsg <- 32611

# define the file you want to work with
# this is the hyperspectral flightline from the hard drive
f <- paste0(dataDir, "/NIS1_20130612_104651_atmcor.h5")

# define clip.extents
clip.extent <- flight1.plots

# calculate extent of H5 file
h5.ext <- create_extent(f)
h5.ext

## class       : Extent 
## xmin        : 296670 
## xmax        : 297409 
## ymin        : 4098356 
## ymax        : 4103734

# turn the H5 extent into a polygon to check overlap
h5.ext.poly <- as(extent(h5.ext), 
                  "SpatialPolygons")

crs(h5.ext.poly) <- crs(clip.extent)

# test to see that the extents overlap
gIntersects(h5.ext.poly, 
            clip.extent)

## [1] TRUE

# Use the clip extent to create the index extent that can be used to slice out data from the 
# H5 file
# xmin.index, xmax.index, ymin.index, ymax.index
# all units will be rounded which means the pixel must occupy a majority (.5 or greater)
# within the clipping extent

index.bounds <- vector("list", length(clip.extent))

index.bounds <- calculate_index_extent(extent(clip.extent),
                                h5.ext)
index.bounds

## [1]  359  462 2815 3656

# this is what i wrote to a csv!